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A modified Prckopa's approach is considered for the problem of optimum allocation in 
multivariate stratified random sampling. An example is solved by applying the proposed 
methodology. 

Key Words: Multivariate stratified random sampling, stochastic programming, opti- 
' mum allocation, integer programming, chance constraints. 

o ' 

Mathematics Subject Classification: 62D05, 90Cf5, 90C29, 90CfO 

O ' 1 Introduction 

One of the statistical tools most commonly used in many fields of scientific research is the 
theory of probabilistic sampling. In diverse practical situations, the probabilistic model 
of stratified random sampling is frequently applied. Although there are different ways to 
allocate the sample in strata, the optimum allocation has been found to be a useful approach, 
see (Stuart, 1954), Cochran (1977), Sukhatme et al. (1984) and Thompson (1997). 

From a multivariate point of view, there are, basically, two approaches for solving the 
problem of optimum allocation in multivariate stratified random sampling. When a cost 
function is denned as the objective function subject to certain functions of variances to 
be within a given region, the problem of the optimum allocation in multivariate strati- 
fied random sampling is stated as a deterministic uniobjective mathematical programming 
problem, see Arthanari and Dodge (1981) among others. Alternatively, when the objective 
function is defined as some functions of variances subject to cost restrictions, the problem 
has been proposed implicitly and explicitly as a deterministic multiobjective mathematical 
programming problem, see Cochran (1977), Sukhatme et al. (1984) and Dfaz-Garcfa and 
Ulloa (2008). 

On the other hand, Prekopa (1978) considers the approach wherein population variances 
are random variables and formulated the corresponding optimum allocation problem as a 
stochastic (or probabilistic) mathematical programming problem, termed specifically chance 
constraints approach, see Charnes and Cooper (1963). Namely, Prekopa (1978) minimizes 
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a cost function subject to inequality restrictions in terms of the estimated variance of the 
stratified mean of each characteristic, restrictions that are allowed to be violated with 
certain probability. An alternative approach suggested by Diaz-Garcia and Ulloa (2008) is 
developed by Kozak and Wang (2010) from a stochastic point of view. 

This work states the optimum allocation in multivariate stratified random sampling 
as a stochastic integer programming problem, specifically, a modified Prekopa's approach 
is proposed. Section 2 includes some notation and definitions on multivariate stratified 
random sampling and summarizes properties on the asymptotic normality of the sample 
covariance matrices. The optimum allocation in multivariate stratified random sampling 
via chance constraints methodology is studied in Section 3. Finally an application of the 
approach is presented in Section 4. 



2 Preliminaries on multivariate stratified random sampling 

Consider a population of size N, divided into H sub-populations (strata). We wish to find 
a representative sample of size n and an optimum allocation rule for the strata, meeting 
the following requirements: i) to minimize the variance of the estimated mean, subject to 
a budgetary constraint; or ii) to minimize the cost subject to a constraint on the variances; 
this is the classical problem in optimum allocation in univariate stratified sampling, see 
Cochran (1977), Sukhatme et al. (1984) and Thompson (1997). However, if more than 
one characteristic (variable) is being considered, then the problem is known as optimum 
allocation in multivariate stratified sampling. For a formal expression of the problem of 
optimum allocation in multivariate stratified sampling, consider the following notation. 

2.1 Notation 

The subindex h = 1, 2, . . . , H denotes the stratum, i = 1,2,..., Nh or the unit within 
stratum h and j = 1,2, ... ,G denotes the characteristic (variable). Moreover: 

iV/j Total number of units within stratum h 

rih Number of units from the sample in stratum h 



Y h = (Yl...Y%) 

= (Y/ji . . . Y/jjv h 

yh = (yi ■ ■ ■ Yh) 

= (yhi ■ ■ -yhn h )' 



Nh x G matrix population in stratum h; Yhi is the 
i-th G-dimensional value of the i-th unit in stratum h 
rih x G matrix sample in stratum h; yhi is the 



i-th element of the G-dimensional random sample 
in stratum h 

y J hi Value obtained for the i-th unit in stratum h 

of the j-th characteristic 



(m, . . . , uh)' Vector of the number of units in the sample 

Nh 
N 



Wh = -rr Relative size of stratum h 



-j 1 Nh 

Yh = ~TT~^~]yhi Population mean in stratum h of the j-th characteristic 



N h - 

(Y h , . . . , Y^)' Population mean vector in stratum h 



rl ttG, 
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1 nh 



y h = (vl, ■ ■ ■ >V%Y 



K 



E^ 



y ST (y 
s„ 



— i 

ST ' 



Coy(y ST ) 
Cov(y ST ) 



Sample mean in stratum h of the j-th characteristic 
Sample mean vector in stratum h 

Estimator of the population mean in multivariate 

stratified sampling for the j-th characteristic 
Estimator of the population mean vector in 
multivariate stratified sampling 
Variance-covariance matrix in stratum h 



1 — — 

Sh = T^E( y '" _ Y h)i.Yhi - Y h )' 



N h = 

where Sh - k is the covariance in stratum h of the 
j-th and fc-th characteristics; furthermore 
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j-22(Vhi-yh)(vhi-Vh), and 



i h 



Vh) 



Estimator of the covariance matrix in stratum 



Sh 



rih 



1 

— E( yw -y h )(yhi-y h Y 



where Sh jk is the sample covariance in stratum h of the 
j-th and fc-th characteristics; furthermore 



Sh, 



rih 



— ^2(yhi~Vh)(yhi -V$), and 



„2 



i=l 
1 



■n-h 



1 = 1 



(vL - vi) 



Variance-covariance matrix of y 
Estimator of the covariance matrix of y , 

J ST ' 

it is denoted as Cov(y ST ) = Cov(y ), and defined as 

1 ST 



Cov(yl T ,t ST ) Var(/ ST ) ••• Cov(^ T ,»f T ) 



V Cov(y« T ,yi T ) Cov(y° T ,y 2 ST ) 



Var(y« T ) / 



= E 



W h 2 s h 



E 



W h Sh 
N 



Estimated covariance of and y where 

cSv(t ST ,y J ST ) = Cov(iJ^ T ), with 



Cov(y J ST ,y 



ST 1 



h=l 

Var(y 3 ST 



in, 



h = l 



N 



ii 

E 



Wh 2 s 2 hj 
n h 



E 



W h s 



hj 



N 



h=l "'" h=l 

Cost per G-dimensional sampling unit in stratum h and let 

c = (ci, . . .,c G y. 



Where if a £ !R G , a' denotes the transpose of a. 
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2.2 Asymptotic normality 

Now, the asymptotic distribution of the estimator, s^, of the covariance matrix is stated. 
First, consider the following notation and definitions. 

A detailed discussion of operator "vec" , "vech" , Moore- Penrose inverse, Kronecker prod- 
uct, commutation matrix and duplication matrix may be found in Magnus and Neudecker 
(1988), among many others. For convenience, some notation is introduced, although in 
general it adheres to standard notation. 

For all matrix A there exists a unique matrix A + which is termed Moore-Penrose inverse 
of A. 

Let A be an m x n matrix and Bapxg matrix. The mp x nq matrix defined by 



aiiB 



011B 



011B 



is termed the Kronecker product (also termed tensor product or direct product) of A and 
B and written A B. Let C be an m x n matrix and Cj its j-th column, then vec C is the 
mn x 1 vector 

vec C 



The vectors vec C and vec C clearly contain the same mn components, but in different 
order. Therefore there exists a unique mn x mn permutation matrix which transforms 
vec C into vec C This matrix is termed the commutation matrix and is denoted K mn (If 
m = n, it is often written K n instead of K mn ). Hence 

K mn vec C = vec C'. 

Similarly, let B be a square nxn matrix. Then vechB (also denoted as v(B)) shall denote 
the n(n + l)/2 x 1 vector that is obtained from vecB by eliminating all supradiagonal 
elements of B. If B = B', vechB contains only the distinct elements of B, then there exists 
a unique n 2 x n(n + l)/2 matrix termed duplication matrix, which is denoted by D n , such 
that D n vechB = vecB and D+ vecB = vechB. Finally, note that (vechB)' = vech'B. 

Now, with the above mathematical tools and based in the extension given in Hajek 
(1961), the multivariate version of Hajek's theorem is restated in terms of sampling theory 
terminology, which is explained in detail in Diaz Garcia and Ramos-Quiroga (2011). 



Lemma 2.1. 



Let S u be a G x G symmetric random matrix defined as 



£(y* - Y„)(y* - Y u )'. 



Suppose that for A = (Ai, . . . , \k)' , any vector of constants, k = G(G + l)/2, 



X'(Mt -vech S y vech' A > e max A^e^ (M^ — vech S y vech' S„) 

1 < a < k L 



(i) 
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where = (0, . . . , 0, 1, 0, . . . , 0)' is the a-th vector of the canonical base of e > and 
independent of v > 1 and 



G 



J2(y»i - Y v )(y wi - Y v )' ® (y* - Y u )(y vi - Y,)' 



i=l 



D 



+' 
G ' 



is the fourth central moment. Assume that n u — > oo, N u — n v — > oo, N u — > oo, and that, 
for all j = 1,...,G, 



lim [ — — 



> lim 



max > 

l<n<---<«n I/ <A r ^ J 



7/ J _T^Y_Q 2 



(2) 



where 



rn 



i=l 



Then, vechHj, is asymptotically normal distributed as 

vechEi, — > A4(E(vech S^), Cov(vech S v )), 

with 



E(vech S v ) 



n,, 



n v -\ 



vech S„, 



and 



Covfvech S u ) = - (M^ — vech S„ vech' S^) . 



(3) 
(4) 



K - 1) ; 

n v is the sample size for a simple random sample from the v-th population of size N u . 
Then: 

Theorem 2.1. Under assumptions in Lemma 2.1, the sequence of sample covariance ma- 
trices s u are such that vechSj, has an asymptotical normal with asymptotic mean and co- 
variance matrix given by (3) and (4), respectively. 

Proof. This follows immediately from Lemma 2.1, only observe that 



n, 



— j- ^{Yui - Y v ){Yvi - y v )' 



i=l 
n v 

n v — \ 



(y„-Y„)(y„-Y„)', 



where 



n,, - 1 



1 and (y„ — Y u )(y u — Y u )' — > in probability. 



□ 



□ 
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Remark 2.1. Observe that it is possible to found the asymptotic distribution of vecS„, 
but this asymptotic normal distribution is singular, because Cov(vec S„) is singular. This is 
due to the fact Cov(vec S^) is the G 2 x G 2 covariance matrix in the asymptotic distribution 
of vec S u and, because is symmetric, then vec Sj, has repeated elements. In this case, 
vec S,, is asymptotically normally distributed as (see Muirhead (1982)) 

vecSj, -4- A/" G 2(E(vecS^), Cov(vec E„)), 



where 



E(vec 3 U ) 



n,. 



n v — \ 



vec S u , 



Cov(vec S t 



n. 



in v - 1)5 



(9Jt^ — vec Sj/ vec' S u ) , 



and 



9JT = — 

v N u 



~ Y -)(y- - Y -)' ® (y* - Y -)(y- - Y -)' 

.i=l 



□ 



The following assertion is an immediate consequence of Theorem 2.1. 
Theorem 2.2. Let Cov(y ST ) be the estimator of the covariance matrix ofy ST , then 

vechCov(y ST ) = - ^) vechs ft 

h=l 

is asymptotically normally distributed; furthermore 

vech Cov(y ST ) -4 A4 ( E (vech Cov(y ST ) J , Cov (vech Cov(y s , 



where 



E (vech Cov(y 



ST 



h=l 



N J n h - 1 



vechS^, 



(5) 
(6) 



Cov (vechCov(y ST ) 



(7) 



and 



= — D+ 



N h 



J2&hi ~ Y ft )(y w - Y h )' ® (y M - Y fc )(yM - Y ft )' 



i=i 



D 



+' 
G ■ 



Finally, note that the asymptotically normal distributions of vech S^, vecCov(y ST ) and 
vechCov(y ST ) are in terms of the parameters Y^, vechS/j, 9JT^ and M^; then, from Rao 
(1973, iv), pp. 388-389), approximations of asymptotic distributions can be obtained, mak- 
ing the following substitutions 



h -» Yhi 



vech Sh — > vech s/j,, 9Jt^ — > rti/j and — > m^ 



(8) 



6 



where 



and 



5^(y« - yh)(y^ - yj' ® (ym - yj(y« - yJ' 



™ 4 



i=l 



D 



G ' 



J^(yw - Yh)(yhi - y h Y ® (yw - y/»)(yiw - yj' 



i=l 



3 Modified Prekopa's approach 

Optimum allocation in multivariate stratified random sampling was proposed as the follow- 
ing deterministic mathematical programming problem 



min c n + Co 

n 

subject to 
y^(yi T )<4, j = l,2,...,G 
2<n h <N h , h = l,2,...,H 

n h G N, 



(9) 



where Vq are desired precisions assigned to the variances of the sample mean Var(y^ r ), 
j = 1,2, ... ,G. This approach has been treated in detail by Arthanari and Dodge (1981). 

From a stochastic point of view of (9), Prekopa (1978) proposes the following chance 
constraints mathematical program 



min c'n + cq 

n 

subject to 

P (Var(^ ST ) < v j ) > P0 , j = l,2,...,G 
2<n h <N h , h = l,2,...,H 

n h E N, 



(10) 



where < po < 1 is a specified probability. 

The present work considers the following alternative chance constraints mathematical 
programming problem 



min c n + cq 

n 

subject to 

P (Cov(y ST ) < A) > po 

2<n h <N h , h = l,2,...,H 
n h e N 



(11) 



vech Cov(y ) -4 Af k ( E ( vech Cov(y ) ) , Cov ( vech Cov(y c 



where A > is a constant matrix. 

From Dfaz-Garcfa and Ulloa (2008), note that Cov(y ST ) is an explicit function of n, 
and so it must be denoted as as Cov(y sT ) = Cov(y ST (n)). In addition, assume that 
Cov(y ST (n)) is a positive definite matrix for all n, Cov(y ST (n)) > 0. Now, let ni and 
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ri2 be two possible values of the vector n and, recall that, for A and B positive definite 
matrices, A>B<^A — B>0. Hence, there exists a function / such that: / : S — > 9ft, 

Cov(y ST (m)) < Cov(y ST (n 2 )) & f (^(y^m))) < / (Cov(y ST (n 2 ))) (12) 

with Cov(y ST (n)) S S C ^ G ^ G+l )l 2 and S is the set of positive definite matrices. 

Then, (11) can be reduced to the following chance constraints mathematical program 

min c'n + cq 

n 

subject to 

p (/ (cov(y ST )) < r) > Po 

2<n h <N h , h = l,2,...,H 
n h e N 

vechCov(y ) 4 Af k (E (vechCov(y ) ) , Cov ( vechCov(y 



(13) 



There are many possibilities for the definition of /(•), see Dfaz-Garcfa and Ulloa (2008). In 
particular, it is of interest when / = tr (Cov(y ST )^ . Among many others options, it is also 

interesting the case when / = Cov(y sT ) in (13) which is described in detail in Section 4, 
although its application in a real problem poses some algorithmic and numerical challenges 
still under study. 



4 Application 

Lets consider the results of a forest survey conducted in Humboldt County, California, 
originally reported in Arvanitis and Afonja (1971). The population was subdivided into 
nine strata on the basis of the timber volume per unit area, as determined from aerial 
photographs. The two variables included in this example are the basal area (BA) 1 in 
square feet, and the net volume in cubic feet (Vol.), both expressed on a per acre basis. 
The variances, covariances and the number of units within stratum h are listed in Table 1. 
For this example, the matrix optimisation problem under approach (13) is 

min n'c + cq 

n 

subject to 

V{c°^yl T ) VarGSJ )- )- P0 

(14) 

2^n h = 1000 

h=l 

2<n h <N h , h = l,...,9 
Cov(y ST ) 4 Mxa (E (Cov(y ST )) ,Cov (vecCcv(y ST ))) 

n h e N. 

1 In forestry terminology, 'Basal area' is the area of a plant perpendicular to the longitudinal axis of a 
tree at 4.5 feet above ground. 
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Table 1: Variances, covariances and the number of units within each stratum 



Variance 



Stratum 


N h 


„ a 
Ch 


BA 


Vol. 


Covariance 


1 


11 131 


2.5 


1 557 


554 830 


28 980 


2 


65 857 


3.0 


3 575 


1 430 600 


61 591 


3 


106 936 


1.5 


3 163 


1 997 100 


72 369 


4 


72 872 


2.5 


6 095 


5 587 900 


166 120 


5 


78 260 


2.0 


10 470 


10 603 000 


293 960 


6 


51 401 


2.0 


8 406 


15 828 000 


357 300 


7 


24 050 


2.5 


20 115 


26 643 000 


663 300 


8 


46 113 


3.0 


9 718 


13 603 000 


346 810 


9 


102 985 


3.5 


2 478 


1 061 800 


39 872 



"These are simulated costs, also Co is taken as 



4.1 Solution when /(•) = tr( ) 
Observe that by (5), (6) and (7) 

tr Cov (y ST ) ~ TV (E (tr Cov (y ST )) , Var (tr Cov {y ST ))) 

where 



G H 



E trCov(y ST ) =££ 



W h 2 W h \ n h 



,r^\ n h N J n h -l 
j=i h=i 



o2 



G H 



V, ( b Cov(y OT ) ) = g g - « - (Si?) . 



and 



N h 



T70 
h 



i=l 



Standardising the function / in equation (14), it is seen that 



tr Cov (y gT ) - E(tr Cov (y ST )) r - E(tr Cov (y 5T )) 



with 



Var(trCov(y 5T )) y Var(tr Cov (y ST )) 

r - E(tr Cov (j ST )) 



> Po, 



Po = $ 



Var(tr Cov (ygy)) t 



where $(•), denotes the standard normal distribution function. Let e po be the value of the 
standard normal random variable such that 3>(e Po ) = po, in such way that the inequality 
can be established as 

r-E(trCov(y 5r _))^ . 



<h 



Var(tr Cov (y 5T )) 



Ysr))) 
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which holds only if 

r - E(trCov (y ST )) 



— e po> 



_'Var(trCov(y 5T )) 

or equivalently 



E(tr Cov (y ST )) + e po y Var(tr Cov (y ST )) - r < 0. (15) 

Hence, taking into account (8), the equivalent deterministic problem to the stochastic math- 
ematical programming (14), is given by 

min n'c + Co 

n 

subject to 



E(tr Cov (y ST )) + e Po yj Var(tr Cov (y ST )) - r < 
9 

= 1000 

h=l 

2<n h <N h , h = l,...,9 
n h 6 N. 



where 



and 



4 1 

»"% = — 

3 n h 



i=l 



( y hi - A 



(18) 



Remark 4.1. Observe that the estimators yi, s\. and m^. of Y h , Sf , and M^. could 
initially be obtained as 

i) results from a pilot (preliminary) sample or 

ii) using the corresponding values of the estimators from another set of variables, X's, 

correlated to the variables Y's. 

It is important to have this in mind in the minimisation step, because for example, the nh 
that appears in expression (18), is the value of (fixed) used in the pilot study. Same 
comment for the expression of the estimator y° h and s\.. While the n^'s that appear in 
expressions (16) and (17) are the decision variables. □ 
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4.2 Solution when /(•) = | • | 

Assume the following alternative stochastic matrix mathematical programming problem 

min n'c + cq 



subject to 
P (/ (Co^(y ST )) <r) > Po 

9 

n h = 1000 

h=l 

2<n h <N h , h = l,...,9 



vechCov(y ST ) ->■ M G xG (vechO GxG , Cov (vech Cov(y si / 

n h G N, 



where Cov(y ) 



vech 



vcchCov(y ST ) - E (vcchCov(y ST ) 



and vech 1 is the inverse function of function vech. 
Then, the restriction in (19), is 



Cov(y ST ) 



< r) >p 



Cov I vech Cov(y t 



B <g> B = N, 



which for G = 2 and assuming that Cov ( vech Cov (y qT ) ) is such that 

implies that 
where 



Cov(y c 



< tINI 1 /^ > 



Po 



tvt ^fW h 2 W h y n h . 4 , > 



n h N J (n h -l) 



™ 4 



i=l 



see Remark 4.1, and 



vI/frlNl 1 / 4 ), 



with denotes the distribution function of Cov(y ST ) , see Delannay and Caer (2000). 

Let r po be the percentile of a random variable such that VP(r Po ) = po, in such way that the 
inequality can be established as 



(19) 



which holds only if 



Cov(y__) < TIN) 1 / 4 ) > ¥( 



INI 1 / 4 > 



'Poll 



PO ' 
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where the density of Z = Cov(y ST ) is, see Delannay and Caer (2000) 
dG{z) 1 



9 z {z) = -y= exp(z) 1 - erf {V2z 



dz » zK 1 ^2 
where erf(-) is the usual error function defined as 



z > 0, 



2 f x 

erf(x) = —i= / exp(— t 2 )dt. 
V 71 " Jo 



Thus, by (8), the equivalent deterministic problem to the stochastic mathematical program- 
ming problem (19), is given by 

min n'c + Co 

n 

subject to 

9 



1000 



h=l 

2<n h <N h , h = l,...,9 
n h G N, 

Table 2 includes the optimum allocation for each characteristic, BA and Vol (the second 
and third rows) from a deterministic point of view. Also appear (on fourth and fifth rows) 
the optimal allocations via the deterministic problem (9), identified in the table with the 
name Prekopa, and the deterministic version of (13) when /(•) = tr(-). These results are 
presented in their stochastic version in the 7-10th rows. The last three columns show the 
minimum values of the individual variances for the respective optimum allocations and the 
cost identified by each method. The results were computed using the commercial software 
Hyper LINGO/PC, release 6.0, see Winston (1995). The default optimisation methods used 
by LINGO to solve the nonlinear integer optimisation programs are Generalised Reduced 
Gradient (GRG) and branch-and-bound methods, see Bazaraa et al. (2006). Finally, note 
that, for this sampling study, there is not a great discrepancy between the different methods 
among the sizes of the strata. And for the multivariate solutions, the biggest cost difference 
appears in the deterministic version of Prekopa's method. 

5 Conclusions 

There is a vast literature on the problem of sample allocation in stratified sampling. A 
natural approach considers a cost minimisation problem subject to variance restrictions. 
This paper follows Prekopa's approach by setting the problem into the area of stochastic 
optimization. It is recognized that this is a more realistic approach because, in general, the 
population variances of the strata are unknown and therefore requires estimating them. As 
a result, problem (9) really falls within the scope of stochastic mathematical programming 
which incorporates the inherent uncertainty of estimators in a natural way. 

The approach is not without its drawbacks, it is difficult to give general rules for electing 
the value function /(•), potentially there are an infinite number of possibilities. In this paper 
we have chosen to work with /(A) = |A| and /(A) = tr(A) which can be interpreted as 
a generalised variance and as an average variance respectively. However, the responsibility 
for the selection or definition of that function, lies wholly with the expert in the field of 
application. 
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Table 2: Sample sizes and estimator of variances for the different allocation rules 



Allocation 


ni 


ni 


U3 


714 


715 


ng 


717 


ns 


ng 




Vai(j£ T ) 


Cost 










Deterministic approach 










BA° 


10 


78 


171 


123 


194 


114 


75 


90 


91 


5.599 


5766.161 


2225.5 


Vol 6 


6 


51 


139 


123 


204 


163 


90 


109 


61 


6.502 


5499.996 


2194.0 


Prckopa c 


10 


78 


171 


123 


194 


111 


75 


90 


91 


5.599 


5766.161 


2225.5 


trCo^(y ST ) d 


6 


17 


127 


114 


186 


149 


80 


102 


59 


7.071 


5992.921 


2014.0 








Stochastic approach^ 












BA 


10 


79 


168 


125 


196 


117 


76 


91 


95 


5.939 


5693.354 


2248.0 


Vol 


6 


48 


129 


113 


189 


150 


82 


102 


60 


6.988 


5933.759 


2034.0 


Prckopa 


11 


79 


169 


123 


196 


117 


78 


91 


96 


5.921 


5680.571 


2034.0 


trCov^J 


6 


48 


129 


114 


188 


151 


81 


102 


60 


6.988 


5933.752 


2034.0 



"With i)i = 6 
6 With vl = 6000 

c With vl = 6 and With u<j = 6000 
d With t = 6000 
e With po = 0.50 
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